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Abstract 

In this paper, we investigate synchronization of coupled second-order linear har- 
monic oscillators with random noises and time delays. The interaction topology is 
modeled by a weighted directed graph and the weights are perturbed by white noise. 
On the basis of stability theory of stochastic differential delay equations, algebraic 
graph theory and matrix theory, we show that the coupled harmonic oscillators can 
be synchronized almost surely with perturbation and time delays. Numerical examples 
are presented to illustrate our theoretical results. 
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1. Introduction 

Synchronization, as an emergent collective phenomenon of a population of units with 
oscillatory behaviors, is one of the most intriguing in nature and plays a significant role 
in a variety of disciplines such as biology, sociology, physics, chemistry and technology 
[21 [151 [21]. O ne celebrated model for synchronization is the Kuramoto model [9], which 
is mathematically tractable and described by a system of structured ordinary differential 
equations. The original Kuramoto formulation assumes full connectivity of the network, 
i.e. whose interaction topology is a complete graph. Recent works generalize the Kuramoto 
model to nearest neighbor interaction and the underlying topologies may be general com- 
plex networks, see e.g. [H ITU] ITS]. 

Another classical model for synchronization is the harmonic oscillator network [3] ITTl 
123] . which is the subject of this paper. Recently, Ren [T7] investigates synchronization 
of coupled second-order linear harmonic oscillators with local interaction. Due to the 
linear structure, the ultimate trajectories to which each oscillator converges over directed 
fixed networks are shown explicitly and milder convergence conditions than those in the 
case of Kuramoto model can be derived. Consensus problems, known as the ability of 
an ensemble of dynamic agents to reach a common value in an asymptotical way, are 
prominent applications of synchronization in engineering and computer science; see [H] 
for a recent survey. The harmonic oscillator networks are related to the second-order 
consensus protocols addressed in [18] EOj [22], [2l] . Compared with these works, where the 
consensus equilibrium for the velocities of agents is a constant, the positions and velocities 
are synchronized to achieve oscillating motion by utilizing harmonic oscillator schemes 
(c.f. Remark 4 below). 

Since noise is ubiquitous in nature, technology, and society, the motion of oscillator is 
inevitably subject to disturbance in the environments. On the other hand, time delay is 
also unavoidable because the information spreading through a complex network is charac- 
terized by the finite speeds of signal transmission over a distance. Although random noise 
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and time delay have been studied extensively in exploring synchronization (or consensus 
problems) by means of theoretical and numerical methods, they have seldom been analyt- 
ically treated in synchronization of coupled harmonic oscillators. Motivating this idea, the 
objective of this paper is to deal with synchronization conditions for coupled harmonic os- 
cillators over general directed topologies with noise perturbation and communication time 
delays. The main tools used here are algebraic graph theory matrix theory and stochastic 
differential delay equation theory. 

The rest of the paper is organized as follows. In Section 2, we provide some preliminar- 
ies and present the coupled harmonic oscillator network model. In Section 3, we analyze 
the synchronization stability of this model and give sufficient conditions for almost surely 
convergence. Numerical examples are given in Section 4 to validate our theoretical results. 
Finally, the conclusion is drawn in Section 5. 

2. Problem formulation 

By convention, R represents the real number set; /„ is an n x n identity matrix. For 
any vector x, x T denotes its transpose and its norm ||x|| is the Euclidean norm. For a 
matrix A, denote by \\A\\ the operator norm of A, i.e. \\A\\ = sup{||^4x|| : ||x|| = 1}. Re(z) 
denotes the real part of z G C. 

Throughout the paper we will be using the following concepts on graph theory (see 
e.g. |5j) to capture the topology of the network interactions. 

Let Q = (V, £ , .A) be a weighted directed graph with the set of vertices V = {1, 2, • • • , n} 
and the set of arcs f C VxV. The vertex i'm.Q represents the ith oscillator, and a directed 
edge (i, j) G £ means that oscillator j can directly receive information from oscillator i. 
The set of neighbors of vertex i is denoted by Mi = {j G V| (j, i) G £}. A = (oy) G M. nxn 
is called the weighted adjacency matrix of Q with nonnegative elements and > if 
and only if j G M%. The in-degree of vertex i is defined as d{ = S?=i a ij- The Laplacian 
of Q is defined as L = D — A, where D = diag(di,d2, ■ •• ,d n ). A directed graph Q is 
called strongly connected if there is a directed path from i to j between any two distinct 
vertices i, j G V. There exists a directed path from vertex i to vertex j, then j is said to 
be reachable from %. For every vertex in directed graph Q, if there is a path from vertex 
i to it, then we say i is globally reachable in Q. In this case, we also say that Q has a 
directed spanning tree with root i. 

Consider n coupled harmonic oscillators connected by dampers and each attached to 
fixed supports by identical springs with spring constant k. The resultant dynamical system 
can be described as 



where x% G M. denotes the position of the ith oscillator, k serves as a positive gain, and ay- 
characterizes interaction between oscillators i and j as mentioned before. 

Here we study a leader-follower version (c.f. Remark 1) of the above system, and more- 
over, communication time delay and stochastic noises during the propagation of informa- 
tion from oscillator to oscillator are introduced. In particular, we consider the dynamical 
system of the form: 
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Xi(t) + kxiit) + 22 a ij(^i(* - r ) - %j(t - T )) + h{xi(t - t) - x (t - r)) 

+ [ a ^(^i( t ~ T ) ~ ^J'C* ~ T )) + Pi(^t(* - T) - 2Co(< - T))]lVi(t) = 0, i = 1, ■ ■ • ,n, 

(2) 

scb(t) + fc»o(t) =0, (3) 

where r is the time delay and xq is the position of the virtual leader, labeled as oscillator 0, 
which follows Equation (|3j) describing an undamped harmonic oscillator. We thus concern 
another directed graph Q D Q associated with the system consisting of n oscillators and 
one leader. Let B = diag(&i,--- ,b n ) be a diagonal matrix with nonnegative diagonal 
elements and 6« > if and only if G Mi- Let W(t) := (wi(t), ■ ■ ■ ,w n (t)) T be an n- 
dimensional standard Brownian motion. Hence, Wi{t) is one-dimensional white noise. To 
highlight the presence of noise, it is natural to assume that <7« > if j G Mi, and o"jj = 
otherwise; pij > if j G A/i, and pij = otherwise. Also let A a = (o"y) G M nxn and 
S CT = diag(pi, • • • , p n ) be two matrices representing the intensity of noise. Moreover, let 
a i = YTj=\ °ih D ° = diag(cii, • • • , u n ), and L a = D a - A a . 

Remark 1. Consensus problems of self-organized groups with leaders have broad appli- 
cations in swarms, formation control and robotic systems, etc.; see e.g. [7| 1 13^ . In 
multi-agent systems, the leaders have influence on the followers' behaviors but usually in- 
dependent of their followers. One therefore transfers the control of a whole system to that 
of a single agent, which saves energy and simplifies network control design Q \B$- Most of 
the existing relevant literatures assume a constant state leader, while our model serves to 
be an example of oscillating state leader on this stage. 

Let ri = Xi and Vi = x"i for i = 0, 1, • • • , n. By denoting r = (ri,--- ,r n ) T and 
v = (vi, ■ ■ ■ , v n ) T , we can rewrite the system ([2]), ([3]) in a compact form as: 

dr(i) = v(t)dt, (4) 
dv(t) = [-kr(t)-(L + B)v(t-T) + Bv (t-T)l]dt 

+ [-(L a + B a )v(t - t) + B a v (t - r)l] dW, (5) 

dr (f) = v (t)dt, dv (t) = -kr (t)dt, (6) 

where 1 denotes an n x 1 column vector of all ones (with some ambiguity; however, the 
right meaning would be clear in the context). 

Remark 2. Note that Vi depends on the information from its in-neighbors and itself. In 
the special case that time delay r = and A a = B a = 0, algorithms {^])-([6jj are equivalent 
to algorithms (12) and (13) in \17j . 

3. Convergence analysis 

In this section, the convergence analysis of systems ([l])-([6]) is given and we show that 
n coupled harmonic oscillators (followers) are synchronized to the oscillating behavior of 
the virtual leader with probability one. 

Before proceeding, we introduce an exponential stability result for the following n- 
dimensional stochastic differential delay equation (for more details, see e.g. j6j) 
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dx(t) = \Ex{t) + Fx(t - r)]dt + g(t, x(t),x(t - r))dW(t), (7) 

where E and F are n x n matrices, g : [0, oo) x W 1 x R n — > R nxm which is locally 
Lipschitz continuous and satisfies the linear growth condition with g(t,0, 0) = 0, VF(t) is 
an m-dimensional standard Brownian motion. 

Lemma l.QllJ) Assume that there exists a pair of symmetric positive definite n x n 
matrices P and Q such that P(E + F) + (E + F) T P = —Q. Assume also that there exist 
non-negative constants a and (3 such that 

tr&ce[g T (t, x, y)g(t, x, y)} < a\\x\\ 2 + (3\\y\\ 2 (8) 

for all (t,x,y) G [0, oo) x R n x M n . Let \ m i n (Q) be the smallest eigenvalue of Q. If 



(a + 0)\\P\\ + 2\\PF\\ V2t(4t(||£|| 2 + ||F|| 2 ) + a + J) < A min (Q), 

then the trivial solution of Equation fiTty is almost surely exponentially stable. 

We need the following lemma for Laplacian matrix. 

Lemma 2.Q19J) Let L be the Laplacian matrix associated with a directed graph Q . Then 
L has a simple zero eigenvalue and all its other eigenvalues have positive real parts if 
and only if Q has a directed spanning tree. In addition, LI = and there exists p S W 1 
satisfying p > 0, p T L = and p T l = 1. 

Let 

r (t) := cos(Vkt)r (0) + \ sin(V^>o(0), 
v (t) := -Vksm(Vkt)r (0) + cos(Vkt)v (0) . 

Then it is easy to see that ro(t) and vo(t) solve ©. Let r* = r — rol, v* = v — Vol. 
Invoking Lemma 2, we can obtain an error dynamics of (jU)-© as follows 

de(t) = [Ee(t) + Fe(t - r)]dt + He(t - r)dW(t), (9) 

where 

^ E=( ° T J « V F=(" „V H ° 



V* J ' V ~ kI n / \ -L-B J ' V -At - B, 

and W(t) is an 2n-dimensional standard Brownian motion. 
Now we present our main result as follows. 

Theorem 1. Suppose that vertex is globally reachable in Q . If 



\\Hf\\P\\ + 2||PF|| y/8r 2 [(k V l) 2 + ||F|| 2 ] + 2r||#|| 2 < A min (Q), (10) 

where k V 1 := max{/c, 1}, P and Q are two symmetric positive definite matrices such that 
P{E + F) + (E + F) T P = -Q, then by using algorithms (^)-([2|), we have 

r{t) - r (t)l -» 0, v(t) - vo(t)l -» 

almost surely, as t — > oo. Here, ro and Wo are given as above. 

Proof. Clearly, it suffices to prove the trivial solution e(t; 0) = of © is almost surely 
exponential stable. 
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Let {Xi : i = 1, • • • , n} be the set of eigenvalues of — L — B. Since vertex is globally 
reachable in Q, from Lemma 2 it follows that —L — B is a stable matrix, that is, Re(Aj) < 
for all i. 

Let n be an eigenvalue of matrix E+F and tp = (ipf , if2) T be an associated eigenvector. 
We thus have 

In \ / <Pl \ ( <Pl 



-kl n -L - B ) \ tp 2 J V (f2 



which yields (— L — B)ipi = ^ ^ k ( fi and tp% / 0. Hence \x satisfies fi 2 — Aj/i + k = 0. The 



2n eigenvalues of E + F are shown to be given by fii± = v 2 for i = 1, ■ ■ ■ , n. Since 

Re(Aj) < 0, we get Re(/Xj_) = Re(^ — ~ ) < for z = 1, • • • , n. From = k 

it follows that Hi + and are symmetric with respect to the real axis in the complex 
plane. Accordingly, Re( / Uj+) < for i = 1, • • • ,n; furthermore, + F is a stable matrix. 
By Lyapunov theorem, for all symmetric positive definite matrix Q there exists a unique 
symmetric positive definite matrix P such that 

P(E + F) + (E + F) T P = -Q. (11) 

On the other hand, we have trace {e T H T He) < ||-£f || 2 ||e|| 2 . Therefore, ([8]) holds with 
a = and f3 = \\H || 2 . Note that \\E\\ = k V 1. We then complete our proof by employing 
Lemma 1. □ 

Remark 3. Note that the result of Theorem 1 is dependent of the choice of matrices P 
and Q. From computational points of view, the solution to Lyapunov matrix equation Ul\) 
may be expressed by using Kronecker product; \\H\\ = \\L a + B a \\ and \\F\\ = \\L + B\\ hold. 

Remark 4. The algorithms |^P-(0|) can also be applied to synchronized motion coordi- 
nation of multi- agent systems, as indicated in \1T$ Section 5. 

When deviations between oscillator states exist, we may exploit the following algorithm 
to take the place of Equation ©: 

dv(t) = [-k(r(t)-5)-(L + B)v(t-T) + Bv {t-T)l]dt 

+ [-(L a + B a )v(t - r) + B a v (t - r)l]dW, (12) 

where 5 = (5\, ■ ■ ■ , 5 n ) T is a constant vector denoting the deviations. Similarly, we obtain 
the following result. 

Corollary 1. Suppose that vertex is globally reachable in Q, and condition U0\) holds, 
then by using algorithms Q), (0) and IWi) . we have 

r(t) - S - r Q (i)l -► 0, v(t) - v (t)l -»■ 

almost surely, as t — > oo. Here, r^ and vq are defined as in Theorem 1. 

4. Numerical examples 

In this section, we provide numerical simulations to illustrate our results. 

We consider a network Q consisting of five coupled harmonic oscillators including one 
leader indexed by and four followers as shown in Fig. 1. We assume that a$j = 1 if 
j E Mi and aij = otherwise; b% = 1 if £ Mi and bi = otherwise. Note that vertex 
is globally reachable in Q. For simplicity, we take the noise intensity matrices L a = 0.1L 
and B a = 0.1B. We take Q = 1% with \ m i n (Q) = 1. By straightforward calculation, it 
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is obtained that \\H\\ = 0.2466 and \\F\\ = 2.4656. Two different gains k are explored as 
follows: 

Firstly, we take k = 0.6 such that ||J5|| = 1 > k. We solve P from Equation (|lip and 
get ||P|| = 8.0944 and ||-P-F|| = 4.1688. Hence the condition (|10p in Theorem 1 is satisfied 
by taking time delay r = 0.002. Thus, the oscillator states are synchronized successfully as 
shown in Fig. 2 and Fig. 3 with initial values given by e(0) = (—5, 1, 4, —3, —8, 2, —1.5, 3) T . 

Secondly, we take k = 2 such that ||J5|| = k > 1. In this case we obtain ||P|| = 8.3720, 
||P-F|| = 7.5996 and the condition (|10p is satisfied by taking time delay r = 0.001. Thereby 
the oscillator states are synchronized successfully as shown in Fig. 4 and Fig. 5 with the 
same initial values given as above. 

We see that the value of k not only has an effect on the magnitude and frequency of 
the synchronized states (as implied in Theorem 1), but also affects the shapes of synchro- 
nization error curves ||r*|| and ||f*||. 

5. Conclusion 

This paper is concerned with synchronization of coupled harmonic oscillators with 
stochastic perturbation and time delays. Based on the stability theory of stochastic dif- 
ferential delay equations, we have shown that the coupled second-order linear harmonic 
oscillators are synchronized (i.e. follow the leader) with probability one provided the 
leader is globally reachable and the time delay is sufficient small. Numerical simulations 
are presented to illustrate our theoretical results. Since we only investigate the case when 
the network topology is fixed, how to consider the time varying topology is our future 
research. 
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Fig. 1 Directed network Q for five coupled harmonic oscillators involving one leader. 
Q has — 1 weights. 

Fig. 2 Synchronization error ||r*|| for k = 0.6 and r = 0.002. 
Fig. 3 Synchronization error ||u*|| for k = 0.6 and r = 0.002. 
Fig. 4 Synchronization error ||r*|| for k = 2 and r = 0.001. 
Fig. 5 Synchronization error ||i>*|| for k = 2 and r = 0.001. 
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